## Endorsements from Republican politicians can increase confidence in U.S. elections 
## Katherine Clayton and Robb Willer
## Last updated: December 12, 2022


# Initial settings --------------------------------------------------------

library(tidyverse)
library(estimatr)
library(scales)
library(mediation)
library(broom)
library(magrittr)
library(TOSTER)
rm(list = ls())


# Read data ---------------------------------------------------------------

df <- read.csv("data/clayton_willer.csv")

# PCA ---------------------------------------------------------------------

legit <- df %>% 
  dplyr::select(legit20_nation, legit20_rightful, legit20_result) %>% 
  na.omit()

legit.pc <- prcomp(legit)
summary(legit.pc)
ltm::cronbach.alpha(legit)

trust <- df %>% 
  dplyr::select(trust_scale, trust_nation, trust_fairly) %>% 
  na.omit()

trust.pc <- prcomp(trust)
summary(trust.pc)
ltm::cronbach.alpha(trust)

strength.pid <- df %>% 
  dplyr::select(partysoc1_pre, partysoc2_pre) %>% 
  na.omit()

strength.pid.pc <- prcomp(strength.pid)
summary(strength.pid.pc)
ltm::cronbach.alpha(strength.pid)

strength.pid.post <- df %>% 
  dplyr::select(partysoc1_post, partysoc2_post) %>% 
  na.omit()

strength.pid.post.pc <- prcomp(strength.pid.post)
summary(strength.pid.post.pc)
ltm::cronbach.alpha(strength.pid.post)

rep.favor <- df %>% 
  dplyr::select(favor_rep, trust_rep) %>% 
  na.omit()

rep.favor.pc <- prcomp(rep.favor)
summary(rep.favor.pc)
ltm::cronbach.alpha(rep.favor)


# Build scales and rescale variables --------------------------------------

df.cleaned <- df %>% 
  rowwise() %>% 
  mutate(legit.out = mean(c(legit20_nation, legit20_rightful, legit20_result), na.rm = TRUE),
         trust.out = mean(c(trust_scale, trust_nation, trust_fairly), na.rm = TRUE),
         pid.soc.mod = mean(c(partysoc1_pre, partysoc2_pre), na.rm = TRUE),
         pid.soc.out = mean(c(partysoc1_post, partysoc2_post), na.rm = TRUE),
         rep.favor.out = mean(c(favor_rep, trust_rep), na.rm = TRUE)) %>% 
  ungroup() %>% 
  mutate(trump.out = rescale(favor_trump_post, to = c(0, 1)),
         legit.out = rescale(legit.out, to = c(0, 1)),
         trust.out = rescale(trust.out, to = c(0, 1)),
         pid.soc.out = rescale(pid.soc.out, to = c(0, 1)),
         rep.favor.out = rescale(rep.favor.out, to = c(0, 1))) %>% 
  mutate(pid.soc.mod2 = NA,
         pid.soc.mod2 = ifelse(pid.soc.mod >= median(pid.soc.mod), 1, pid.soc.mod2),
         pid.soc.mod2 = ifelse(pid.soc.mod < median(pid.soc.mod), 0, pid.soc.mod2)) %>% 
  mutate(trump.mod = NA,
         trump.mod = ifelse(favor_trump_pre >= median(favor_trump_pre), 1, trump.mod),
         trump.mod = ifelse(favor_trump_pre < median(favor_trump_pre), 0, trump.mod)) %>% 
  mutate(pid7 = NA,
         pid7 = ifelse(pid == 1 & pid7_rep == 1, 1, pid7),
         pid7 = ifelse(pid == 1 & pid7_rep == 2, 2, pid7),
         pid7 = ifelse(pid == 3 | pid == 4 & pid7_ind == 1, 3, pid7)) %>% 
  mutate(male = ifelse(gender == 1, 1, 0),
         white = ifelse(race == 1, 1, 0),
         college = ifelse(education >= 5, 1, 0))

# Sample details ----------------------------------------------------------

mean(df.cleaned$Duration..in.seconds.)
median(df.cleaned$Duration..in.seconds.)
prop.table(table(df.cleaned$qual_watch))
prop.table(table(df.cleaned$qual_sound))
prop.table(table(df.cleaned$humorous))
prop.table(table(df.cleaned$white))
prop.table(table(df.cleaned$gender))
prop.table(table(df.cleaned$age))
prop.table(table(df.cleaned$college))
prop.table(table(df.cleaned$income))


# Pooling test ------------------------------------------------------------

control <- df.cleaned %>% filter(group == "control")
placebo <- df.cleaned %>% filter(group == "placebo")

t.test(control$legit.out, placebo$legit.out)
t.test(control$trust.out, placebo$trust.out)


# Combine placebo and control ---------------------------------------------

df.cleaned <- df.cleaned %>% 
  mutate(group = ifelse(group == "placebo", "control", group)) %>% 
  mutate(group = factor(group, levels = c("control", "dem", "rep")))


# Main hypothesis tests and visualization ---------------------------------

h1 <- estimatr::lm_robust(legit.out~group + age + male + white + college + income, df.cleaned)
summary(h1)
summary(multcomp::glht(h1, linfct = c("grouprep - groupdem = 0"))) 

h2 <- estimatr::lm_robust(trust.out~group + age + male + white + college + income, df.cleaned)
summary(h2)
summary(multcomp::glht(h2, linfct = c("grouprep - groupdem = 0"))) 

h1.tidy <- h1 %>%
  broom::tidy() %>% 
  filter(term == "groupdem" | term == "grouprep") %>% 
  mutate(sig = ifelse(p.value < 0.05, TRUE, FALSE)) %>% 
  mutate(outcome = "legitimacy")

h2.tidy <- h2 %>%
  broom::tidy() %>% 
  filter(term == "groupdem" | term == "grouprep") %>% 
  mutate(sig = ifelse(p.value < 0.05, TRUE, FALSE)) %>% 
  mutate(outcome = "trust")

out1 <- rbind(h1.tidy, h2.tidy)

out1$group <- factor(out1$outcome, levels = c("legitimacy","trust"), labels = c("Perceived 2020 legitimacy", "Trust in US elections"))

dv_lab <- c(rep("Perceived 2020 legitimacy", 2), rep("Trust in US elections", 2))
mod <- c(rep("h1", 2), rep("h2", 2))
null.value <- c(rep(0, 4))

out1 <- cbind(out1, dv_lab, mod, null.value) %>% 
  mutate(contrast = NA,
         contrast = ifelse(term == "groupdem", "Democrat cue", contrast),
         contrast = ifelse(term == "grouprep", "Republican cue", contrast))

df.viz <- out1 %>% 
  mutate(estimate = as.numeric(estimate)) %>% 
  mutate(lab = ifelse(abs(estimate) < .01, round(estimate, digits = 3), round(estimate, digits = 2))) %>% 
  mutate(lab = as.character(lab)) %>% 
  mutate(lab = str_replace(lab, fixed("0."), ".")) %>% 
  mutate(lab = str_c(lab, 
                     c("***", "**", "*", "") %>% 
                       magrittr::extract(
                         p.value %>% 
                           findInterval(
                             c(-Inf, .001, .01, .05, Inf)
                           )
                       ) )) %>% 
  mutate(dv_lab = factor(dv_lab, c("Trust in US elections", "Perceived 2020 legitimacy"))) %>% 
  mutate(std.error = as.numeric(std.error),
         conf.low = as.numeric(conf.low),
         conf.high = as.numeric(conf.high)) %>% 
  mutate(contrast = factor(contrast, levels = c("Republican cue", "Democrat cue")))

main.viz <- ggplot(data = df.viz) +
  geom_vline(
    aes(xintercept = 0),
    linetype = "dashed") + 
  geom_linerange(
    aes(xmin = conf.low, 
        xmax = conf.high, 
        y = dv_lab, 
        color = contrast,
        group = contrast),
    data = df.viz,
    position = position_dodge2(width = .4), 
    show.legend = F
  ) +
  geom_point(
    aes(estimate, dv_lab,
        color = contrast,
        group = contrast),
    data = df.viz,
    position = position_dodge2(width = .4),
    shape = 21,
    size = 15,
    fill = "white"
  ) +
  geom_point(
    aes(estimate, dv_lab,
        group = contrast),
    data = df.viz,
    position = position_dodge2(width = .4),
    shape = 21,
    size = 7.75,
    color = "white",
    fill = "white", 
    show.legend = F
  ) +
  xlim(-.1, .1) + 
  geom_text(
    aes(
      estimate, dv_lab, label = lab,
      color = contrast,
      group = contrast
    ),
    data = df.viz,
    position = position_dodge2(width = .4),
    size = 4,
    show.legend = F
  ) +
  scale_y_discrete(
    breaks = c(
      "Trust in US elections",
      "Perceived 2020 legitimacy"
    ),
    labels = c(
      "Trust in US elections",
      "Perceived 2020 legitimacy"
    ) %>%
      str_wrap(width = 20)
  ) +
  scale_color_manual(
    values = c("#2fabe1", "#e91b23") %>% rev,
    guide = guide_legend(
      override.aes = list(shape = 16, size = 3),
      reverse = T
    )
  ) +
  facet_grid(
    contrast~., 
    scales = "free_y",
    space = "free_y",
    labeller = label_wrap_gen(width = 15)
  ) +
  labs(
    x = "Difference relative to control",
    y = "") +
  theme(
    legend.position = "none",
    panel.background  = 
      element_rect(color = "grey95",
                   fill = "grey95"),
    panel.grid = element_blank(),
    strip.background = element_rect(color = "white",
                                    fill = "white"),
    legend.background = element_rect(color = "white",
                                     fill = "white"),
    panel.grid.minor = element_blank(),
    axis.text=element_text(size=12),
    axis.title=element_text(size=14),
    strip.text.y = element_text(size = 12)
  ) 

main.viz

ggsave("figures/main_viz.png", main.viz, width = 10, height = 6)


# Exploratory analysis: Effects on unscaled outcomes ----------------------

legit20_nation.mod <- estimatr::lm_robust(legit20_nation~group + age + male + white + college + income, df.cleaned)
summary(legit20_nation.mod)
summary(multcomp::glht(legit20_nation.mod, linfct = c("grouprep - groupdem = 0"))) 

legit20_rightful.mod <- estimatr::lm_robust(legit20_rightful~group + age + male + white + college + income, df.cleaned)
summary(legit20_rightful.mod)
summary(multcomp::glht(legit20_rightful.mod, linfct = c("grouprep - groupdem = 0"))) 

legit20_result.mod <- estimatr::lm_robust(legit20_result~group + age + male + white + college + income, df.cleaned)
summary(legit20_result.mod)
summary(multcomp::glht(legit20_result.mod, linfct = c("grouprep - groupdem = 0"))) 

trust_scale.mod <- estimatr::lm_robust(trust_scale~group + age + male + white + college + income, df.cleaned)
summary(trust_scale.mod)
summary(multcomp::glht(trust_scale.mod, linfct = c("grouprep - groupdem = 0"))) 

trust_nation.mod <- estimatr::lm_robust(trust_nation~group + age + male + white + college + income, df.cleaned)
summary(trust_nation.mod)
summary(multcomp::glht(trust_nation.mod, linfct = c("grouprep - groupdem = 0"))) 

trust_fairly.mod <- estimatr::lm_robust(trust_fairly~group + age + male + white + college + income, df.cleaned)
summary(trust_fairly.mod)
summary(multcomp::glht(trust_fairly.mod, linfct = c("grouprep - groupdem = 0"))) 


# Causal mediation analysis -----------------------------------------------

df.med <- df.cleaned %>%
  filter(!(group == "dem")) %>% 
  filter(!is.na(group)) %>% 
  mutate(treat = ifelse(group == "rep", 1, 0))

med.fit.robust <- estimatr::lm_robust(manip_1~treat + age + male + white + college + income, df.med)
summary(med.fit.robust)

med.fit <- lm(manip_1~treat + age + male + white + college + income, df.med)
out.fit1 <- lm(legit.out~manip_1 + treat + age + male + white + college + income, df.med)
out.fit2 <- lm(trust.out~manip_1 + treat + age + male + white + college + income, df.med)

set.seed(05055)

med.out1 <- mediate(med.fit, out.fit1, treat = "treat", mediator = "manip_1", robustSE = TRUE, sims = 1000)
summary(med.out1)

med.out2 <- mediate(med.fit, out.fit2, treat = "treat", mediator = "manip_1", robustSE = TRUE, sims = 1000)
summary(med.out2)

medsens(med.out1, rho.by = 0.1, effect.type = "indirect") %>% summary()
medsens(med.out2, rho.by = 0.1, effect.type = "indirect") %>% summary()

# Moderators --------------------------------------------------------------

summary(estimatr::lm_robust(legit.out~group + trump.mod + group:trump.mod + age + male + white + college + income, df.cleaned))
summary(estimatr::lm_robust(trust.out~group + trump.mod + group:trump.mod + age + male + white + college + income, df.cleaned))

summary(estimatr::lm_robust(legit.out~group + pid.soc.mod2 + group:pid.soc.mod2 + age + male + white + college + income, df.cleaned))
summary(estimatr::lm_robust(trust.out~group + pid.soc.mod2 + group:pid.soc.mod2 + age + male + white + college + income, df.cleaned))

summary(estimatr::lm_robust(legit.out~group + factor(pid7) + group:factor(pid7) + age + male + white + college + income, df.cleaned))
summary(estimatr::lm_robust(trust.out~group + factor(pid7) + group:factor(pid7) + age + male + white + college + income, df.cleaned))


# Additional exploratory outcomes -----------------------------------------

rep.favor <- lm_robust(rep.favor.out~group + age + male + white + college + income, df.cleaned)
summary(rep.favor)
summary(multcomp::glht(rep.favor, linfct = c("grouprep - groupdem = 0"))) 

trump.out <- lm_robust(trump.out~group + trump.mod + age + male + white + college + income, df.cleaned)
summary(trump.out)
summary(multcomp::glht(trump.out, linfct = c("grouprep - groupdem = 0"))) 

pid.soc.out <- lm_robust(pid.soc.out~group + pid.soc.mod2 + age + male + white + college + income, df.cleaned)
summary(pid.soc.out)
summary(multcomp::glht(pid.soc.out, linfct = c("grouprep - groupdem = 0"))) 


# Equivalence tests -------------------------------------------------------

dm.repfavor <- df.cleaned %>% 
  group_by(group) %>%   
  summarise(mean = mean(rep.favor.out, na.rm = TRUE), 
            sd = sd(rep.favor.out, na.rm = TRUE), 
            n = n(),
            se = sd(rep.favor.out, na.rm = TRUE)/sqrt(n())) %>% 
  mutate(lower_ci = mean - qt(1 - ((1 - 0.95) / 2), n - 1) * se) %>% 
  mutate(upper_ci = mean + qt(1 - ((1 - 0.95) / 2), n - 1) * se)

TOSTtwo.raw(m1 = dm.repfavor$mean[3],
            m2 = dm.repfavor$mean[2],
            sd1 = dm.repfavor$sd[3],
            sd2 = dm.repfavor$sd[2],
            n1 = dm.repfavor$n[3],
            n2 = dm.repfavor$n[2],
            low_eqbound = -0.03,
            high_eqbound = 0.03,
            alpha = 0.05, var.equal = FALSE)


TOSTtwo.raw(m1 = dm.repfavor$mean[3],
            m2 = dm.repfavor$mean[1],
            sd1 = dm.repfavor$sd[3],
            sd2 = dm.repfavor$sd[1],
            n1 = dm.repfavor$n[3],
            n2 = dm.repfavor$n[1],
            low_eqbound = -0.03,
            high_eqbound = 0.03,
            alpha = 0.05, var.equal = FALSE)


# Additional analyses -----------------------------------------------------

table(control$favor_biden)
